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The ground state phase diagram of the two-dimensional attractive Hubbard model with popula- 
tion imbalance is explored using a mean field ansatz. A linear programming algorithm is used to 
identify the blocked states, such that the population imbalance can be imposed exactly. This allows 
to explore regions of the number-projected phase diagram that can not be obtained with the conven- 
tional Bogoliubov-de Gennes ansatz. The Fulde-Ferrel-Larkin-Ovchinnikov (FFLO) phase of pairs 
with non-zero momentum is found to be the ground state over a wide range of parameters, while 
phase separation occurs only in a limited region at small population imbalance. Through a particle- 
hole transformation these results can be related to the underdoped repulsive Hubbard model, where 
the FFLO phase takes the form of a particle-hole condensate that exhibits a spontaneous restoration 
of spin symmetry. 
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The Bardeen, Cooper and SchriefFer (BCS) theory of 
superconductivity [1] has been very successful in explain- 
ing not only the properties of superconductors but also 
many properties of atomic nuclei and in recent years 
also properties of ultracold Fermi gases. In this latter re- 
spect a particularly interesting setting is offered by opti- 
cal lattices, because they allow an unprecedented control 
over the interaction parameters and geometry of strongly 
interacting quantum systems [3]. The most important 
degrees of freedom of such a system are well described 
by the minimal Hubbard Hamiltonian f3] , 

^ a+flj,, + J7^a+a+a,xa,|. (1) 

The first summation runs over nearest neighbors only. 
The Hubbard interaction parameter U reflects the short- 
range interaction between particles at the same lattice 
site. In the atomic physics setting a^^ creates an atom 
in the optical lattice site i and the spin tr distinguishes 
between two different types of atoms or between two dis- 
tinct hyperfine levels [5J . In condensed matter physics 
creates an electron with spin-projection a in the Wannier 
orbital centered at lattice site i. This model is relevant 
for high-Tc superconductors where the superconductivity 
is believed to originate from electrons or holes moving in 
two-dimensional cupper-oxide lattices [6J. 

For a system with an equal number of spin-up and 
spin-down particles and a negative U , i.e. an attractive 
interaction, particles of opposite spin pair up to form 
zero-momentum pairs that can condense |7]. If there is 
an asymmetry in the particle numbers of both spins, e.g. 
an excess of spin-up particles, then not every spin-up par- 
ticle will find a spin-down partner to pair up. It is not 
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yet fully understood what happens in that case. On the 
one hand the unpaired particles want to reduce their ki- 
netic energy, but on the other hand they block part of the 
phase space where the system wants to develop pairing 
correlations. These conflicting tendencies have inspired 
many scenarios [8]: the simple BCS picture, where the ex- 
cess particles minimize their kinetic energy but block the 
pairing correlations up to some extent; breached pairing, 
where the excess partices are pushed to higher kinetic en- 
ergies in order to permit full pairing correlations at the 
Fermi surface of the minority component or where some 
particles of the minority component move to higher ki- 
netic energies in order to develop pairing correlations at 
the Fermi surface [51 [TD]; or deformed Fermi surfaces, 
such that minority and majority Fermi surfaces coincide 
in certain directions [TT]. Already in the sixties Fulde 
and Ferrel and independently Larkin and Ovchinnikov 
suggested another possibility: non-zero momentum pairs 
such that pairing correlations can develop between pairs 
at different kinetic energies, and thus bridge the gap be- 
tween the Fermi surfaces of the two species [12j [T3j [14] . 
Even though convincing arguments have been made that 
such a state might be energetically favorable [THl [H] , it 
has been elusive to demonstrate experimentally, and only 
recently experiments have provided some evidence for its 
existence in the quasi-twodimensional heavy fermion su- 
perconductor CeCoIns [17 . Another possibility is that 
the excess particles form unpaired regions in real space 
[181 119] . The first experiments with unbalanced mix- 
tures of ultracold fermions in traps exhibited such phase 
separation [50} [H]. Which of the above scenarios pre- 
vails will depend sensitively on the population imbalance, 
the interaction strength and the shape of the kinetic en- 
ergy dispersion curve. Particularly, it has been shown 
that the FFLO scenario is the preferred configuration in 
the strongly interacting limit [22]. Here I present a new 
algorithm to determine the blocked states at the mean 
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field level. This method is then applied to the attrac- 
tive Hubbard model in two dimensions. Previous studies 
have already indicated that in the Hubbard model with 
imbalanced spin populations the FFLO phase gains sta- 
bility near half-filling because non-zero momentum pairs 
can maximally exploit the enhanced phase space ofl^ered 
by the nesting of the Fermi surfaces and the correspond- 
ing Van Hove singularity in the middle of the band |23| . 
The ground state phase diagram in this regime is of par- 
ticular importance because it might be accessible in ex- 
periments with ultracold atoms in optical lattices and 
it can be related to the phase diagram of the repulsive 
Hubbard model through a unitary particle-hole transfor- 
mation, which might shed light on the pairing mechanism 
in High-Tc superconductors 

Mean-field equations can be obtained from a varia- 
tional ansatz with a trial state \'^) for which Wick's 
theorem holds, such that the expectation value of 
the interaction can be split into a sum of den- 
sity X density terms U{^\a'j'^ai^\^){^\af^aii\^) - 
U{^\af^aii\^){9\af^a^ij\^) (the Hartree terms) and pair- 
ing terms U{^\dl-^al,\^){^\aiiai^\^). A variation on the 
trial state leads to the Bogohubov-de Gennes equations, 
also known as the Hartree- Fock-Bogoliubov equations. A 
restriction of the formalism to states that do not mix 
spin-up and spin-down components and to the Fulde- 
Ferrel modulation, (^'|a+^a+Jv]/) = e^'i'-'Aq/C/, is still 
general enough to compare BCS, breached and FFLO 
phases and phase separation and to identify the relevant 
degrees of freedom. Minimization of the energy to \^) 
for iV| spin-up and spin-down particles in a L x i 
lattice, leads to a mean field Hamiltonian of the form 

Hnif,q = ^ [(Ck + rfk) Ck|CkT + (ek - rfk) Cq_k|Cq_ki 



+A,c+ c+_k^ A;cq_kiCkT +9N^N^ + ^ (2) 
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with g = -U/L'^, ek - (ek + Cq-k - 5(^^T + ^i)) /2 " M 
and dk (ck — eq_k + g{N^ — A^i)) /2. The eigenprob- 
lem of Eq. ^ for the unblocked momentum states reduces 
to the BCS form, except that the pairs now have a finite 
momentum q. The eigenstates are solutions of the BCS- 
like equations 



k unblocked 



1 



(3) 



where 
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"k.T = "-q-k,i 

menta, and or 
Lagrange multiplier 



-f A^. The occupation numbers are 
(1 ~ ek/ii'k) /2 for the unblocked mo- 
1 for the blocked momenta. The 
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to be adjusted such that 
X^k ("-k.T + '^k,i) = -^T + -^i- imbalance in spin pop- 
ulations, e.g. > Ni, leads to complications due to the 
blocking by unpaired particles [25] . Rather than using a 
different fi„ for each spin and fix A| — Aj^ on average, here 
the difference between the particle numbers is imposed 
exactly by explicitly selecting a set of blocked momenta. 



As shown below this allows to reach certain configura- 
tions that can not be obtained with spin-dependent La- 
grange multipliers. How to select the blocked momenta is 
a subtle question because each configuration of unpaired 
particles leads to another eigenstate of the mean field 
Hamiltonian. For a given value of /i and Aq the optimal 
configuration can be found efficiently using a variation of 
the Needleman-Wunsch algorithm 261 1 ^^^^ known for its 
application to genetic sequence alignment. To this end 
one imagines a graph where the momentum states are 
ordered on the horizontal axis, and the sum of — 
for all previous momentum states appears on the vertical 
axis. For each momentum k one has three options: a pair 
state with energy e\^ — Ek, an unpaired spin-up particle 
with energy ek + dk, and an unpaired spin-down particle 
with energy ek — rfk, represented respectively by a hor- 
izontal jump, a diagonal jump upwards and a diagonal 
jump downwards to the right. There are an exponential 
number of possible paths that lead to the desired pop- 
ulation imbalance A| — A^^, each with a different total 
energy. The optimal path can be split at any momentum 
k: the segment up to that point on the graph necessarily 
has to be the optimal path to get from the origin to that 
point; therefore one has to keep track only of the opti- 
mal paths for each point, and these can be constructed 
by comparing at most three possible paths that arrive in 
each point of the graph. In this way the ground state for 
a given value of Ag can be obtained in O (3(A^| + Ni)L'^) 
operations. By iterating this procedure for each vector q 
on the reciprocal lattice, for several initial values of A^, 
one can determine the absolute ground state. 

The zero temperature phase diagram of the half-filled 
2D attractive Hubbard model was determined by evaluat- 
ing the ground state for a 25 x 25 lattice for N^+N\^ = , 
all possible asymmetries x — {N^ — Ni)/{N^ + -^i) and 
for all possible values of q. Fig|l^ shows the evolution 
of the energy for U = —6. The system starts out in a 
BCS configuration with zero-momentum pairs (B). At 
X — 0.02 the FFLO configurations with q > start 
to take over: first configurations where q aligns along 
the diagonal (FFLO-d), and then at x = 0.19 configu- 
rations with q parallel to a lattice axis (FFLO-p). At 
large asymmetries Aq goes to zero and everything con- 
verges to the normal state (N). For q > this point is 
shifted to larger asymmetries. This means that FFLO 
mechanism allows the system to avoid the the Clogston- 
Chandrasekhar limit [371 HI]- Instead, the system un- 
dergoes a crossover from the FFLO region to the nor- 
mal state at a larger asymmetry, that depends on the 
interaction strength U. The minimal energy curve in 
Figjlji displays a concave behavior for x < 0.3. A grand 
canonical approach based on a difference in Lagrange 
multipliers, h = (/i| — /i|)/2, would result in first order 
transition at he = 2.077 from the balanced BCS con- 
figuration to the FFLO-p minimum at Xc — 0.28, and 
would never produce a ground state with particle num- 
bers corresponding to an asymmetry x in between and 
Xc- This means that the results in this range can only be 
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FIG. 1: Mean field results for the half filled 25 x 25 attrac- 
tive Hubbard model with population imbalance: (a) Energy 
per site shifted by hcX, for various configurations at f7 = —6; 
(b) number projected ground state phase diagram. Configu- 
rations N, B, PS, FFLO-d and FFLO-p explained in the text. 
The momentum densities for the points marked by x and + 
are shown in Fig.Q. 



obtained by imposing the exact particle number differ- 
ence. The FFLO states in this region could be stabilized 
by long range spatial modulations, as explained already 
by Larkin and Ovchinnikov |13) . and by the residual in- 
teraction. Apart from the N, B, FFLO-d and FFLO-p 
regions, the ground state phase diagram in FigjlJ) also 
shows a region with phase separation (PS), where the 
system splits into balanced BCS domains and fully polar- 
ized domains. The corresponding energy was evaluated 
as = vEf {\N^ - N^\/v) + {I - v)Ep {NJ{1- v)), 
with Efin) the energy of n unpaired spin- up particles in 
the Lx L lattice without paired particles and with Epin) 
the energy of n pairs in the same lattice but without un- 
paired particles, where v was optimized to minimize Epg. 
Note how the value of q evolves as a function of the Hub- 
bard interaction U: q is not a constant function of the 
difference between the two Fermi momenta, as suggested 
by previous mean field studies [24^, but rather depends 
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FIG. 2: Momentum density profiles for the diagonal and par- 
allel FFLO configurations marked in Fig. 0: (a) spm-up par- 
ticles, (b) spin-down particles (c) transformed spin-up holes. 

on U times x as shown by the bands in Fig|T|D. 

Momentum density profiles typical for the FFLO-d and 
FFLO-p regions are shown in Fig|2] They clearly show 
the regions of blocked momenta. Therefore an obvious 
way to confirm the FFLO phases experimentally would 
be to look for such distinctive asymmetries in the spin 
projected momentum densities. Another option would 
be to look for signals in radiofrequency spectroscopy or 
momentum resolved photoemission spectroscopy |29j , al- 
though a more advanced many-body treatment might be 
necessary in order to identify the FFLO phases accu- 
rately. An interesting point about the FFLO phase is 
that the lowest energy degrees of freedom are not given 
by Bogoliubov quasiparticles (they are inhibited by the 
gap Aq), but rather by excitations of the unpaired parti- 
cles, for which there is no gap. The dynamics of these un- 
paired particles is governed by an effective singe-particle 
potential given by dk — E\^ and a p wave interaction of 
first order in U induced by the residual part of the full 
Hamiltonian through the exchange of spin-down particles 
with the pair condensate. The unpaired particles also ex- 
perience an interaction of second order in U , induced by 
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phonons in the pair condensate [3D]. 

Finally, I want to draw attention to the implications of 
these results for the repulsive Hubbard model: one can 
define a particle-hole transformation, Ck| = c^^)_kt' 
Ckx = Ckj, |31j . that converts the half-filled attractive 
model with population imbalance into the underdoped 
repulsive model with equal spin populations. If the above 
calculations give a qualitative picture of the ground state 
for the imbalanced attractive model, then the particle- 
hole transformed equivalent should be relevant for the 
ground state of the underdoped repulsive model. The 
third panel of Fig. ([2| shows the transformed momentum 
profile for the spin-up particles. There is a remarkable 
symmetry with the momentum profile for the spin-down 
particles, that seems to be a general feature of the FFLO 
phase, and that makes sense because there is no a pri- 
ori asymmetry in the repulsive case. The FFLO pair 
condensate transforms into an asymmetric particle-hole 
condensate. The onset of this condensation occurs in 
the same parameter range as the pseudogap |32j in the 
underdoped repulsive model, with an anisotropic gap of 
the order of the Hubbard interaction. The residual in- 
teraction could then produce strong hole-hole correlation 
effects: the anisotropy of the repulsion between the holes 
would enhance the formation of striped patterns, in line 
with the suggestions from Moreo et al. on the relation 
between the attractive FFLO phase and a striped phase 



in the repulsive model |24j . 

Through explicit selection of the blocked momenta, 
the particle number difference A^i — was imposed ex- 
actly in the mean field solution. This has allowed the 
exploration of normal, phase separated, BCS pairing, 
breached pairing and non-zero momentum Fuldc-Fcrrel 
pairing phases in the ground state phase diagram for the 
two-dimensional attractive Hubbard model with popu- 
lation imbalance at half filling. At very small popula- 
tion imbalances, the ground state evolves from normal to 
BCS pairing to phase separation as \U\ is increased. At 
larger imbalances the phase diagram at intermediate in- 
teraction strengths is dominated by non-zero momentum 
FFLO pairs, with a pair momentum q that scales pro- 
portionally to the population imbalance and the interac- 
tion strength. These configurations might be identified 
experimentally from asymmetries in the spin projected 
momentum densities. Through a particle-hole transfor- 
mation these results can be related to the repulsive Hub- 
bard model with equal spin populations below half filling. 
There the FFLO mechanism gives rise to the formation 
of a particle-hole condensate on which the residual inter- 
action builds strong hole-hole correlations. This might 
be related to a pseudogap phase and stripe formation. 

I thank K. M0lmer, N. Nygaard, R. Molina and J. 
Dukelsky for the interesting discussions. 
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